function pca_analysis()
% PCA_ANALYSIS 主成分分析详细实现
% 
% 使用SVD实现PCA的各种变体和应用

fprintf('=== 主成分分析 (PCA) ===\n\n');

%% 1. 基础PCA实现
fprintf('1. 基础PCA实现\n');
fprintf('目标：理解PCA的数学原理\n\n');

% 生成相关数据
rng(42);
n_samples = 200;
n_features = 4;

% 构造具有内在结构的数据
true_components = 2;
W_true = [1, 0.5; 0.8, -0.6; -0.3, 0.9; 0.7, 0.4];  % 真实载荷矩阵
scores_true = randn(true_components, n_samples);
X_clean = W_true * scores_true;

% 添加噪声
noise_level = 0.3;
X_noisy = X_clean + noise_level * randn(n_features, n_samples);

fprintf('数据维度: %d特征 × %d样本\n', n_features, n_samples);
fprintf('真实内在维度: %d\n', true_components);
fprintf('噪声水平: %.1f\n', noise_level);

% 数据中心化
X_mean = mean(X_noisy, 2);
X_centered = X_noisy - X_mean;

fprintf('数据预处理:\n');
fprintf('原始数据均值: [%.3f, %.3f, %.3f, %.3f]\n', X_mean);
fprintf('中心化后均值: [%.2e, %.2e, %.2e, %.2e]\n', mean(X_centered, 2));

% 方法1：协方差矩阵的特征分解
fprintf('\n方法1：协方差矩阵特征分解\n');
C = (X_centered * X_centered') / (n_samples - 1);
fprintf('协方差矩阵:\n');
disp(C);

[V_cov, D_cov] = eig(C);
[eigenvals_cov, idx] = sort(diag(D_cov), 'descend');
V_cov = V_cov(:, idx);

fprintf('特征值: [%.4f, %.4f, %.4f, %.4f]\n', eigenvals_cov);
fprintf('主成分载荷矩阵:\n');
disp(V_cov);

% 方法2：SVD分解
fprintf('\n方法2：SVD分解\n');
[U_svd, S_svd, V_svd] = svd(X_centered, 'econ');
singular_values = diag(S_svd);
eigenvals_svd = (singular_values.^2) / (n_samples - 1);

fprintf('奇异值: [%.4f, %.4f, %.4f, %.4f]\n', singular_values);
fprintf('对应特征值: [%.4f, %.4f, %.4f, %.4f]\n', eigenvals_svd);
fprintf('主成分载荷矩阵:\n');
disp(U_svd);

% 验证两种方法的一致性
fprintf('\n方法一致性验证:\n');
eigenval_diff = norm(eigenvals_cov - eigenvals_svd);
loadings_diff = norm(abs(V_cov) - abs(U_svd), 'fro');  % 符号可能不同
fprintf('特征值差异: %.2e\n', eigenval_diff);
fprintf('载荷矩阵差异: %.2e\n', loadings_diff);

%% 2. 主成分选择
fprintf('\n2. 主成分选择\n');
fprintf('目标：确定保留的主成分数量\n\n');

% 方差解释比
total_variance = sum(eigenvals_svd);
explained_variance_ratio = eigenvals_svd / total_variance;
cumulative_variance = cumsum(explained_variance_ratio);

fprintf('主成分分析结果:\n');
fprintf('成分    特征值    方差解释比    累积解释比\n');
fprintf('----    ------    ----------    ----------\n');
for i = 1:length(eigenvals_svd)
    fprintf('%4d    %6.4f    %10.3f    %10.3f\n', ...
            i, eigenvals_svd(i), explained_variance_ratio(i), cumulative_variance(i));
end

% Kaiser准则（特征值>1）
kaiser_components = sum(eigenvals_svd > 1);
fprintf('\nKaiser准则（特征值>1）: %d个成分\n', kaiser_components);

% 累积方差准则
variance_thresholds = [0.8, 0.9, 0.95, 0.99];
fprintf('累积方差准则:\n');
for threshold = variance_thresholds
    n_components = find(cumulative_variance >= threshold, 1);
    fprintf('保留%.0f%%方差: %d个成分\n', 100*threshold, n_components);
end

% 碎石图分析（寻找"肘部"）
fprintf('碎石图分析:\n');
eigenval_ratios = eigenvals_svd(1:end-1) ./ eigenvals_svd(2:end);
[~, elbow_idx] = max(eigenval_ratios);
fprintf('最大特征值比值在成分%d和%d之间: %.2f\n', ...
        elbow_idx, elbow_idx+1, eigenval_ratios(elbow_idx));

%% 3. 主成分得分和载荷
fprintf('\n3. 主成分得分和载荷\n');
fprintf('目标：分析主成分的含义\n\n');

% 选择前2个主成分
n_components = 2;
loadings = U_svd(:, 1:n_components);
scores = S_svd(1:n_components, 1:n_components) * V_svd(:, 1:n_components)';

fprintf('前%d个主成分分析:\n', n_components);
fprintf('载荷矩阵 (%d×%d):\n', size(loadings));
disp(loadings);

fprintf('得分矩阵统计:\n');
for i = 1:n_components
    fprintf('PC%d: 均值=%.2e, 标准差=%.4f, 范围=[%.2f, %.2f]\n', ...
            i, mean(scores(i,:)), std(scores(i,:)), ...
            min(scores(i,:)), max(scores(i,:)));
end

% 载荷解释
fprintf('\n载荷解释:\n');
feature_names = {'特征1', '特征2', '特征3', '特征4'};
for i = 1:n_components
    fprintf('PC%d的主要贡献特征:\n', i);
    [~, sorted_idx] = sort(abs(loadings(:, i)), 'descend');
    for j = 1:n_features
        fprintf('  %s: %.4f\n', feature_names{sorted_idx(j)}, ...
                loadings(sorted_idx(j), i));
    end
    fprintf('\n');
end

%% 4. 数据重构和压缩
fprintf('\n4. 数据重构和压缩\n');
fprintf('目标：评估不同主成分数量的重构质量\n\n');

components_to_test = 1:n_features;
fprintf('重构质量分析:\n');
fprintf('成分数    重构误差    相对误差    压缩比\n');
fprintf('------    --------    --------    ------\n');

for k = components_to_test
    % k个主成分重构
    X_reconstructed = U_svd(:, 1:k) * S_svd(1:k, 1:k) * V_svd(:, 1:k)' + X_mean;
    
    % 重构误差
    reconstruction_error = norm(X_noisy - X_reconstructed, 'fro');
    relative_error = reconstruction_error / norm(X_noisy, 'fro');
    
    % 压缩比
    original_size = n_features * n_samples;
    compressed_size = k * (n_features + n_samples + 1);
    compression_ratio = compressed_size / original_size;
    
    fprintf('%6d    %8.4f    %8.4f    %6.3f\n', ...
            k, reconstruction_error, relative_error, compression_ratio);
end

%% 5. 异常值检测
fprintf('\n5. 异常值检测\n');
fprintf('目标：使用PCA检测异常样本\n\n');

% 添加一些异常值
n_outliers = 5;
outlier_indices = randperm(n_samples, n_outliers);
X_with_outliers = X_noisy;
X_with_outliers(:, outlier_indices) = X_with_outliers(:, outlier_indices) + ...
                                      3 * randn(n_features, n_outliers);

fprintf('添加%d个异常值\n', n_outliers);

% 重新计算PCA
X_outlier_mean = mean(X_with_outliers, 2);
X_outlier_centered = X_with_outliers - X_outlier_mean;
[U_outlier, S_outlier, V_outlier] = svd(X_outlier_centered, 'econ');

% 使用前2个主成分
k_outlier = 2;
scores_outlier = S_outlier(1:k_outlier, 1:k_outlier) * V_outlier(:, 1:k_outlier)';

% 计算重构误差（SPE统计量）
X_reconstructed_outlier = U_outlier(:, 1:k_outlier) * scores_outlier + X_outlier_mean;
spe_values = sum((X_with_outliers - X_reconstructed_outlier).^2, 1);

% 计算Hotelling T²统计量
eigenvals_outlier = (diag(S_outlier).^2) / (n_samples - 1);
t2_values = sum((scores_outlier ./ sqrt(eigenvals_outlier(1:k_outlier))).^2, 1);

fprintf('异常值检测统计量:\n');
fprintf('SPE阈值 (95%%): %.4f\n', quantile(spe_values, 0.95));
fprintf('T²阈值 (95%%): %.4f\n', quantile(t2_values, 0.95));

% 检测异常值
spe_threshold = quantile(spe_values, 0.95);
t2_threshold = quantile(t2_values, 0.95);

spe_outliers = find(spe_values > spe_threshold);
t2_outliers = find(t2_values > t2_threshold);

fprintf('SPE检测到的异常值: %d个\n', length(spe_outliers));
fprintf('T²检测到的异常值: %d个\n', length(t2_outliers));

% 验证检测效果
spe_detection_rate = length(intersect(spe_outliers, outlier_indices)) / n_outliers;
t2_detection_rate = length(intersect(t2_outliers, outlier_indices)) / n_outliers;

fprintf('检测准确率:\n');
fprintf('SPE方法: %.1f%%\n', 100 * spe_detection_rate);
fprintf('T²方法: %.1f%%\n', 100 * t2_detection_rate);

%% 6. 核PCA
fprintf('\n6. 核PCA (Kernel PCA)\n');
fprintf('目标：处理非线性数据结构\n\n');

% 生成非线性数据
theta = linspace(0, 4*pi, 150);
r = theta / (2*pi);
X_nonlinear = [r .* cos(theta); r .* sin(theta)] + 0.1 * randn(2, 150);

fprintf('非线性数据: %d特征 × %d样本\n', size(X_nonlinear));

% 标准PCA
X_nl_mean = mean(X_nonlinear, 2);
X_nl_centered = X_nonlinear - X_nl_mean;
[U_nl, S_nl, V_nl] = svd(X_nl_centered, 'econ');

fprintf('标准PCA结果:\n');
eigenvals_nl = (diag(S_nl).^2) / (size(X_nonlinear, 2) - 1);
fprintf('特征值: [%.4f, %.4f]\n', eigenvals_nl);
fprintf('方差解释比: [%.3f, %.3f]\n', eigenvals_nl / sum(eigenvals_nl));

% 核PCA（RBF核）
fprintf('\n核PCA (RBF核):\n');
sigma = 1.0;  % 核参数
n_nl = size(X_nonlinear, 2);

% 计算核矩阵
K = zeros(n_nl, n_nl);
for i = 1:n_nl
    for j = 1:n_nl
        K(i, j) = exp(-norm(X_nonlinear(:, i) - X_nonlinear(:, j))^2 / (2 * sigma^2));
    end
end

% 中心化核矩阵
one_n = ones(n_nl, n_nl) / n_nl;
K_centered = K - one_n * K - K * one_n + one_n * K * one_n;

% 特征分解
[V_kernel, D_kernel] = eig(K_centered);
[eigenvals_kernel, idx_kernel] = sort(diag(D_kernel), 'descend');
V_kernel = V_kernel(:, idx_kernel);

% 归一化特征向量
for i = 1:n_nl
    if eigenvals_kernel(i) > 1e-12
        V_kernel(:, i) = V_kernel(:, i) / sqrt(eigenvals_kernel(i));
    end
end

fprintf('核PCA特征值 (前5个): [%.4f, %.4f, %.4f, %.4f, %.4f]\n', ...
        eigenvals_kernel(1:5));

% 核主成分得分
kernel_scores = V_kernel' * sqrt(diag(max(eigenvals_kernel, 0)));

fprintf('核PCA vs 标准PCA:\n');
fprintf('标准PCA第1主成分解释方差: %.3f\n', eigenvals_nl(1) / sum(eigenvals_nl));
fprintf('核PCA第1主成分解释方差: %.3f\n', eigenvals_kernel(1) / sum(eigenvals_kernel));

%% 7. 稀疏PCA
fprintf('\n7. 稀疏PCA\n');
fprintf('目标：获得稀疏的主成分载荷\n\n');

% 生成高维稀疏数据
n_features_sparse = 20;
n_samples_sparse = 100;
true_sparse_components = 3;

% 构造稀疏载荷
W_sparse_true = zeros(n_features_sparse, true_sparse_components);
W_sparse_true(1:5, 1) = [2, 1.5, -1, 0.8, -0.5];      % 第1个成分
W_sparse_true(6:10, 2) = [1.8, -1.2, 0.9, -0.7, 1.1]; % 第2个成分
W_sparse_true(11:15, 3) = [-1.5, 1, -0.8, 1.3, -0.6]; % 第3个成分

% 生成数据
scores_sparse_true = randn(true_sparse_components, n_samples_sparse);
X_sparse = W_sparse_true * scores_sparse_true + 0.2 * randn(n_features_sparse, n_samples_sparse);

fprintf('稀疏PCA数据: %d特征 × %d样本\n', n_features_sparse, n_samples_sparse);
fprintf('真实稀疏成分数: %d\n', true_sparse_components);

% 标准PCA
X_sparse_mean = mean(X_sparse, 2);
X_sparse_centered = X_sparse - X_sparse_mean;
[U_sparse, S_sparse, V_sparse] = svd(X_sparse_centered, 'econ');

fprintf('标准PCA载荷稀疏性:\n');
for i = 1:3
    sparsity_std = sum(abs(U_sparse(:, i)) < 0.1) / n_features_sparse;
    fprintf('PC%d: %.1f%% 接近零的元素\n', i, 100 * sparsity_std);
end

% 简化的稀疏PCA（软阈值）
fprintf('\n稀疏PCA (软阈值方法):\n');
lambda_sparse = 0.3;  % 稀疏性参数

U_sparse_thresholded = U_sparse;
for i = 1:3
    % 软阈值
    U_sparse_thresholded(:, i) = sign(U_sparse(:, i)) .* ...
                                 max(abs(U_sparse(:, i)) - lambda_sparse, 0);
    % 重新归一化
    if norm(U_sparse_thresholded(:, i)) > 0
        U_sparse_thresholded(:, i) = U_sparse_thresholded(:, i) / ...
                                     norm(U_sparse_thresholded(:, i));
    end
end

fprintf('稀疏PCA载荷稀疏性:\n');
for i = 1:3
    sparsity_sparse = sum(abs(U_sparse_thresholded(:, i)) < 1e-6) / n_features_sparse;
    fprintf('PC%d: %.1f%% 零元素\n', i, 100 * sparsity_sparse);
end

% 重构质量比较
X_recon_std = U_sparse(:, 1:3) * (U_sparse(:, 1:3)' * X_sparse_centered) + X_sparse_mean;
X_recon_sparse = U_sparse_thresholded(:, 1:3) * ...
                 (U_sparse_thresholded(:, 1:3)' * X_sparse_centered) + X_sparse_mean;

error_std = norm(X_sparse - X_recon_std, 'fro') / norm(X_sparse, 'fro');
error_sparse = norm(X_sparse - X_recon_sparse, 'fro') / norm(X_sparse, 'fro');

fprintf('重构误差比较:\n');
fprintf('标准PCA: %.4f\n', error_std);
fprintf('稀疏PCA: %.4f\n', error_sparse);
fprintf('稀疏性代价: %.4f\n', error_sparse - error_std);

%% 8. 增量PCA
fprintf('\n8. 增量PCA\n');
fprintf('目标：处理大规模或流数据\n\n');

% 模拟流数据
batch_size = 50;
n_batches = 4;
total_samples = batch_size * n_batches;

fprintf('增量PCA设置:\n');
fprintf('批次大小: %d\n', batch_size);
fprintf('批次数量: %d\n', n_batches);
fprintf('总样本数: %d\n', total_samples);

% 生成所有数据（用于比较）
X_stream_all = randn(n_features, total_samples);
X_stream_all(1:2, :) = 2 * X_stream_all(1:2, :);  % 增强前两个特征的方差

% 批量PCA（基准）
X_all_mean = mean(X_stream_all, 2);
X_all_centered = X_stream_all - X_all_mean;
[U_batch, S_batch, ~] = svd(X_all_centered, 'econ');
eigenvals_batch = (diag(S_batch).^2) / (total_samples - 1);

% 增量PCA
fprintf('\n增量PCA过程:\n');
fprintf('批次    累积样本    前3个特征值\n');
fprintf('----    --------    -----------\n');

% 初始化
X_incremental_mean = zeros(n_features, 1);
C_incremental = zeros(n_features, n_features);
n_seen = 0;

for batch = 1:n_batches
    % 获取当前批次数据
    start_idx = (batch - 1) * batch_size + 1;
    end_idx = batch * batch_size;
    X_batch = X_stream_all(:, start_idx:end_idx);
    
    % 更新统计量
    n_batch = size(X_batch, 2);
    X_batch_mean = mean(X_batch, 2);
    
    if n_seen == 0
        X_incremental_mean = X_batch_mean;
        C_incremental = cov(X_batch');
    else
        % 更新均值
        X_incremental_mean_new = (n_seen * X_incremental_mean + n_batch * X_batch_mean) / ...
                                 (n_seen + n_batch);
        
        % 更新协方差矩阵
        X_batch_centered = X_batch - X_batch_mean;
        C_batch = (X_batch_centered * X_batch_centered') / (n_batch - 1);
        
        delta_mean = X_batch_mean - X_incremental_mean;
        C_incremental = ((n_seen - 1) * C_incremental + (n_batch - 1) * C_batch + ...
                        (n_seen * n_batch / (n_seen + n_batch)) * (delta_mean * delta_mean')) / ...
                        (n_seen + n_batch - 1);
        
        X_incremental_mean = X_incremental_mean_new;
    end
    
    n_seen = n_seen + n_batch;
    
    % 计算当前主成分
    [~, D_incremental] = eig(C_incremental);
    eigenvals_incremental = sort(diag(D_incremental), 'descend');
    
    fprintf('%4d    %8d    [%.4f, %.4f, %.4f]\n', ...
            batch, n_seen, eigenvals_incremental(1:3));
end

% 比较最终结果
fprintf('\n最终结果比较:\n');
fprintf('方法        前3个特征值\n');
fprintf('--------    -----------\n');
fprintf('批量PCA     [%.4f, %.4f, %.4f]\n', eigenvals_batch(1:3));
fprintf('增量PCA     [%.4f, %.4f, %.4f]\n', eigenvals_incremental(1:3));

eigenval_errors = abs(eigenvals_batch(1:3) - eigenvals_incremental(1:3)) ./ eigenvals_batch(1:3);
fprintf('相对误差    [%.2e, %.2e, %.2e]\n', eigenval_errors);

%% 总结
fprintf('\n=== PCA分析总结 ===\n');
fprintf('1. 基础PCA：协方差矩阵特征分解 vs SVD分解\n');
fprintf('2. 主成分选择：Kaiser准则、方差解释比、碎石图\n');
fprintf('3. 载荷和得分：理解主成分的含义\n');
fprintf('4. 数据重构：评估不同成分数的效果\n');
fprintf('5. 异常检测：SPE和Hotelling T²统计量\n');
fprintf('6. 核PCA：处理非线性结构\n');
fprintf('7. 稀疏PCA：获得可解释的载荷\n');
fprintf('8. 增量PCA：处理大规模数据\n\n');

fprintf('PCA应用指南：\n');
fprintf('• 数据预处理：中心化（必须）、标准化（可选）\n');
fprintf('• 成分选择：结合多种准则综合判断\n');
fprintf('• 结果解释：分析载荷矩阵的模式\n');
fprintf('• 异常检测：结合SPE和T²统计量\n');
fprintf('• 非线性数据：考虑核PCA\n');
fprintf('• 高维稀疏：使用稀疏PCA\n');
fprintf('• 大规模数据：采用增量PCA\n\n');

end